↜ Back to index Introduction to Numerical Analysis 2

Solutions of exercises from the previous lecture

See Solutions for Lecture a1

Gaussian elimination (ガウスの消去法)

Our goal in this class is to solve linear systems (1次連立方程式) of the type \begin{aligned} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1,\\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n &= b_2,\\ &\ \ \vdots\\ a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n &= b_n,\\ \end{aligned} where n \in {\mathbb{N}} is the number of equations, a_{11}, a_{12}, …, a_{nn} and b_1, b_2, …, b_n are given numbers and x_1, x_2, …, x_n are the unknowns.

This can be written more compactly as Ax = b, where A is the n\times n matrix \begin{pmatrix} a_{11} & a_{12} &&\cdots & a_{1n}\\ a_{21} & a_{22} & && \vdots\\ a_{31} & a_{32} & a_{33} & \\ \vdots& & & \ddots &\\ a_{n1} & a_{n2} & \cdots & & a_{nn}\\ \end{pmatrix} and x, b \in {\mathbb{R}}^n are the vectors \begin{aligned} b &= (b_1, b_2, \ldots, b_n)\\ x &= (x_1, x_2, \ldots, x_n) \end{aligned}

It is easy to do this if A is upper- or lower triangular.

We therefore need to convert the system Ax=b into a form with a triangular matrix. We can do this rather easily by performing row operations on A and b. Indeed, if we add a multiple of one equation to another equation in the system, we do not change the solution. This is equivalent to adding a multiple of a row of the matrix A to another row, and adding the multiple of the component of b in the same row to the component of b in the other row.

This algorithm is called Gaussian elimination (ガウスの消去法). The following is a basic implementation in Fortran.

implicit none

integer, parameter :: ms = 20    ! maximal allowed n (number of equations)
real a(ms,ms), b(ms), x(ms)      ! arrays for the matrix A, vectors b, x
integer i, j, k, n
real c, s

! read data

! number of equations
read *, n

! read a in column-major order
! the order of loops is important
do j = 1, n
    do i = 1, n
        read *, a(i, j)
    end do
end do

! read b
do i = 1, n
    read *, b(i)
end do

! now we are ready to proceed solving Ax = b

! eliminate columns 1, ..., n-1
do k = 1, n - 1
    if (a(k,k) == 0) then
        print *, "error: a(",k , ", ", k, ") is zero"
        stop
    end if
    do i = k + 1, n
        c = -a(i,k) / a(k,k)
        do j = k, n
            a(i,j) = a(i, j) + c * a(k, j)
        end do
        b(i) = b(i) + c * b(k)
    end do
end do

! add this point time the matrix stored in `a` is a upper triagonal matrix

! compute the solution x of Ax = b
do k = 1, n
    ! we have to go over rows in reverse
    i = n + 1 - k
    s = 0
    do j = i+1,n
        s = s + a(i, j) * x(j)
    end do
    x(i) = (b(i) - s) / a(i, i)
end do

! print the vector x
do i = 1, n
    print *, x(i)
end do

end

You can find versions of this code online on many websites. Our goal is to understand it and develop the skills to write it ourselves based on the explanation of the algorithm and modify it if necessary.

The above code for example fails for a simple matrix A = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}, even though this matrix is invertible and Ax = b has a solution. We will learn how to fix the program later.